library(lmtest)
library(sandwich)
library(foreign)
library(survey)



# Load data ---------------------------------------------------------------

svyz <- read.csv("../data/asylum_data.csv")


# Weight pre-processing ---------------------------------------------------

svyz$weight[svyz$weight > 6] <- 6
svyww <- subset(svyz,!is.na(svyz$weight))


# Country vectors ---------------------------------------------------------

countries<-c("Germany","Hungary","Sweden","Austria","Norway",
             "Switzerland","Denmark","Netherlands","Greece","Czech Republic","Italy",
             "Poland","Spain","France","United Kingdom")

countries2<-c("Germany","Hungary","Sweden","Austria","Norway",
              "Switzerland","Denmark","Netherlands","Greece","Czech Republic","Italy",
              "Poland","Spain","France","United Kingdom","Pooled")



# Percentage who Prefer Prop over Status Quo using Rank Data --------------

svyww$PropOverStatusQuo <- as.numeric(svyww$SharePrank < svyww$ShareFErank)

#First for respondents not assigned to consequences treatment
svywww <- subset(svyww,svyww$ConsTrt == 0)
svywww <- subset(svywww,!is.na(svywww$PropOverStatusQuo))

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est0 <- c()
se0 <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PropOverStatusQuo), data = dfrw, weights = dfrw$weight)
  est0[i] <- svymean(~dfrw.PropOverStatusQuo,surveyobject)[1]
  se0[i] <- sqrt(vcov(svymean(~dfrw.PropOverStatusQuo,surveyobject)))
  
}

est0 <- est0*100
se0 <- se0*100

#Second for respondents assigned to consequences treatment
svywww <- subset(svyww,svyww$ConsTrt == 1)

hold <- svywww
hold$cty <- "Pooled"
svywww <- rbind(svywww,hold)

est1 <- c()
se1 <- c()

for (i in 1:length(countries2)){
  
  dfrw <- subset(svywww,svywww$cty == countries2[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PropOverStatusQuo), data = dfrw, weights = dfrw$weight)
  est1[i] <- svymean(~dfrw.PropOverStatusQuo,surveyobject)[1]
  se1[i] <- sqrt(vcov(svymean(~dfrw.PropOverStatusQuo,surveyobject)))
  
}

est1 <- est1*100
se1 <- se1*100


ddat <- data.frame(countries2=rep(countries2,2),givenshares=rep(c("No","Yes"),each=16),est=c(est0,est1),se=c(se0,se1))

ddat$lower <- ddat$est - 1.96*ddat$se
ddat$upper <- ddat$est + 1.96*ddat$se

ddat


# Delta across Ideo, Pooled across countries, Ideo3 -------------------------------------------------------

svyww$Ideo3 <- NA
svyww$Ideo3[svyww$IdeoScale < 6] <- "1. Left"
svyww$Ideo3[svyww$IdeoScale == 6] <- "2. Center"
svyww$Ideo3[svyww$IdeoScale > 6] <- "3. Right"
ideos <- sort(unique(svyww$Ideo3))


svywww <- subset(svyww,svyww$ConsTrt == 0)

est0 <- c()
se0 <- c()

for (i in 1:length(ideos)){
  
  dfrw <- subset(svywww,svywww$Ideo3 == ideos[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est0A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est0A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est0A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  
  est0B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est0B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)^2)
  
  cov.est0AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)*(dfrw$PreferProp - est0A))
  
  est0[i] <- est0A - est0B
  
  var0 <- vhat.est0A + vhat.est0B - 2*cov.est0AB
  se0[i] <- sqrt(var0)
  
}

est0 <- est0*100
se0 <- se0*100



svywww <- subset(svyww,svyww$ConsTrt == 1)

est1 <- c()
se1 <- c()

for (i in 1:length(ideos)){
  
  dfrw <- subset(svywww,svywww$Ideo3 == ideos[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est1A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est1A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est1A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  est1B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est1B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)^2)
  
  cov.est1AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)*(dfrw$PreferProp - est1A))
  
  est1[i] <- est1A - est1B
  
  var1 <- vhat.est1A + vhat.est1B - 2*cov.est1AB
  se1[i] <- sqrt(var1)
  
}

est1 <- est1*100
se1 <- se1*100


ddat <- data.frame(ideos=rep(ideos,2),givenshares=rep(c("No","Yes"),each=3),est=c(est0,est1),se=c(se0,se1))

ddat$lower <- ddat$est - 1.96*ddat$se
ddat$upper <- ddat$est + 1.96*ddat$se

ddat



# Delta across Want More / Less, Pooled across countries -------------------------------------------------------

svyww$apref <- NA
svyww$apref[svyww$AsylumHome > 0] <- "1. Increase"
svyww$apref[svyww$AsylumHome == 0] <- "2. Same"
svyww$apref[svyww$AsylumHome < 0] <- "3. Decrease"
prefs <- sort(unique(svyww$apref))

svywww <- subset(svyww,svyww$ConsTrt == 0)


est0 <- c()
se0 <- c()

for (i in 1:length(prefs)){
  
  dfrw <- subset(svywww,svywww$apref == prefs[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est0A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est0A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est0A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  
  est0B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est0B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)^2)
  
  cov.est0AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)*(dfrw$PreferProp - est0A))
  
  est0[i] <- est0A - est0B
  
  var0 <- vhat.est0A + vhat.est0B - 2*cov.est0AB
  se0[i] <- sqrt(var0)
  
}

est0 <- est0*100
se0 <- se0*100



svywww <- subset(svyww,svyww$ConsTrt == 1)

est1 <- c()
se1 <- c()

for (i in 1:length(prefs)){
  
  dfrw <- subset(svywww,svywww$apref == prefs[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est1A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est1A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est1A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  est1B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est1B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)^2)
  
  cov.est1AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)*(dfrw$PreferProp - est1A))
  
  est1[i] <- est1A - est1B
  
  var1 <- vhat.est1A + vhat.est1B - 2*cov.est1AB
  se1[i] <- sqrt(var1)
  
}

est1 <- est1*100
se1 <- se1*100


ddat <- data.frame(prefs=rep(prefs,2),givenshares=rep(c("No","Yes"),each=3),est=c(est0,est1),se=c(se0,se1))

ddat$lower <- ddat$est - 1.96*ddat$se
ddat$upper <- ddat$est + 1.96*ddat$se

ddat


# Delta across Knowledge, Pooled across countries -------------------------------------------------------

knows <- sort(unique(svyww$Knowledge))

svywww <- subset(svyww,svyww$ConsTrt == 0)

est0 <- c()
se0 <- c()

for (i in 1:length(knows)){
  
  dfrw <- subset(svywww,svywww$Knowledge == knows[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est0A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est0A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est0A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  
  est0B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est0B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)^2)
  
  cov.est0AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est0B)*(dfrw$PreferProp - est0A))
  
  est0[i] <- est0A - est0B
  
  var0 <- vhat.est0A + vhat.est0B - 2*cov.est0AB
  se0[i] <- sqrt(var0)
  
}

est0 <- est0*100
se0 <- se0*100



svywww <- subset(svyww,svyww$ConsTrt == 1)

est1 <- c()
se1 <- c()

for (i in 1:length(knows)){
  
  dfrw <- subset(svywww,svywww$Knowledge == knows[i])
  
  nw <- nrow(dfrw)
  dfrw$weight <- dfrw$weight/mean(dfrw$weight)
  
  est1A <- (1/nw)*sum(dfrw$PreferProp*dfrw$weight)
  
  dfrw$WeightNorm <- dfrw$weight/sum(dfrw$weight)
  vhat.est1A <- sum(((dfrw$WeightNorm)^2) * (dfrw$PreferProp - est1A)^2)
  
  dfrw$SQPref <- as.numeric(dfrw$SharePref == "1. first entry")
  
  est1B <- (1/nw)*sum(dfrw$SQPref*dfrw$weight)
  vhat.est1B <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)^2)
  
  cov.est1AB <- sum(((dfrw$WeightNorm)^2) * (dfrw$SQPref - est1B)*(dfrw$PreferProp - est1A))
  
  est1[i] <- est1A - est1B
  
  var1 <- vhat.est1A + vhat.est1B - 2*cov.est1AB
  se1[i] <- sqrt(var1)
  
}

est1 <- est1*100
se1 <- se1*100


ddat <- data.frame(knows=rep(knows,2),givenshares=rep(c("No","Yes"),each=4),est=c(est0,est1),se=c(se0,se1))

ddat$lower <- ddat$est - 1.96*ddat$se
ddat$upper <- ddat$est + 1.96*ddat$se

ddat



# Rank Correlation - Support for Prop vs. Win from Prop -------------------


cty.win.order <-c("Germany","Hungary","Sweden","Austria","Norway",
                  "Switzerland","Denmark","Netherlands","Greece","Czech Republic",
                  "Italy","Poland","Spain","France","United Kingdom")
win.order <- seq(1:15)

svywww <- subset(svyww,svyww$InfoTrt == 0 & svyww$ConsTrt == 0)

est.prop <- c()
n <- c()

for (i in 1:length(cty.win.order)){
  
  dfrw <- subset(svywww,svywww$cty == cty.win.order[i])
  
  surveyobject <- svydesign(ids = ~0, variables = data.frame(dfrw$PreferProp), data = dfrw, weights = dfrw$weight)
  est.prop[i] <- svymean(~dfrw.PreferProp,surveyobject)[1]
  
  n[i] <- nrow(dfrw)
  
}


thisdf <- data.frame(cty.win.order,win.order,est.prop,n)

cor.test(thisdf$win.order,thisdf$est.prop,
         alternative = "two.sided",
         method = "spearman")




# Moderation of Information Treatment by Knowledge ------------------------


svywww <- subset(svyww,svyww$ConsTrt == 0)

interaction.est <- c()
interaction.se <- c()

for (i in 1:15){
  out.reg <- lm(PreferProp ~ InfoTrt*Knowledge, data=subset(svywww,svywww$cty == countries[i]), weights = weight)
  print(cty.win.order[i])
  print(coeftest(out.reg, vcov = vcovHC(out.reg, "HC3"))["InfoTrt:Knowledge",])
  interaction.est[i] <- coeftest(out.reg, vcov = vcovHC(out.reg, "HC3"))["InfoTrt:Knowledge","Estimate"]
  interaction.se[i] <- coeftest(out.reg, vcov = vcovHC(out.reg, "HC3"))["InfoTrt:Knowledge","Std. Error"]
}


